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Abstract 

We study, analytically and numerically, mesoscopic fluctuations of the off- 
diagonal matrix elements of the orbital angular momentum between the near- 
est energy levels i = (n x ,n y ) and / = (k x ,k y ) in a rectangular box with 
incommensurate sides. In the semiclassical regime, where the level number 



of i is J\f S> 1, our derivation gives 



L if 



J\f. Numerical simulations, 



using simultaneous ensemble averaging (over the aspect ratios of rectangles) 
and spectral averaging (over the energy interval), are in excellent agreement 
with this analytical prediction. Physically, the mean is dominated by the level 
pairs k x = n x ± 1, k y = n y 1. Also in a rectangular box, we investigate the 
mean orbital susceptibility of a free electron gas and argue that it reduces, 
up to a coefficient, to the two- level van Vleck susceptibility that involves the 
last occupied (Fermi) level i and the first unoccupied level /. This result is 
confirmed numerically as well, albeit the effect of fluctuations is much more 

pronounced for the susceptibility since it is due both to large fluctuations in 

~ 2 \ 

Lij ) and in level separations e j — £i (level bunching) . 
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I. OFF-DIAGONAL MATRIX ELEMENTS OF THE ANGULAR MOMENTUM 



A. Introduction 



Off-diagonal matrix elements of the orbital angular momentum enter into important 
physical quantities, such as magnetic dipole absorption and van Vleck susceptibility. This is 
particularly significant in the situations when the angular momentum is not a good quantum 
number. Such is the case in disordered systems where, in the semiclassical approximation, 
it was shown that in 2D 



Here hp is the wave number of a free particle whose energy ep corresponds to level i, and 
£ and r are, respectively, the mean-free-path and the scattering time due to disorder. This 
result can be derived either by considering the classical magneto-dipole absorption [1], [2] 
or by a direct evaluation [1] using the technique developed in Refs. [3], [4]. 

Disordered systems are classically chaotic. In the semiclassical regime they exhibit "level 
rigidity," which prevents large fluctuations in the level spacings and, in turn, large fluctua- 
tions of the number of levels in an energy interval. [5]- [7] Classically integrable systems, on 
the other hand, exhibit "level bunching" characterized by large fluctuations in level spacings 1 
and high occurrence of small spacings - hence the term - and, in turn, large fluctuations in 
the number of levels in an energy interval. [5] , [6] (It was recently shown that such behavior 
extends only up to a certain energy scale, upon which strong correlations between levels set 
in [8]; Even then, the number level variance exhibit large, non-decaying oscillations around 
the "saturation value.") Consequently, the mesoscopic, or non self-averaging, effects are 
expected to be more pronounced in integrable systems. 



In the absence of resonances, or supersymmetry, such as the case for harmonic oscillator and the 
Kepler problem. 




(1) 
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The fluctuations in the level spacings (and more generally, the specifics of correlations 
between the levels) is, however, just one of the factors contributing to the mesoscopic fluc- 
tuations of physical quantities. Other contributors are expected to fluctuate much stronger 
in classically integrable systems than in classically chaotic systems as well. In the first part 
of this paper we consider the fluctuations of the off-diagonal matrix elements of the angular 
momentum in a rectangular box. We will show, in particular, that 



Li 



if 



~ VJ7 ~ kpL (2) 



where Af >> 1 is the level number of i and L is a rectangle's side. 

Brackets in eq. (1) denote averaging over various realizations of disorder, which is an 
example of "ensemble averaging." A natural extension of the concept of ensemble averaging 
to a rectangular box would be to average over the aspect ratios of the rectangles' sides. 
However, our analytical derivation of eq. (2) is based on "spectral averaging," that is 
averaging over an energy interval. A detailed explanation of our numerical procedure will 
be given in text, but it should be already mentioned that a combined ensemble and spectral 
averaging was performed. The former involves averaging over the aspect ratios chosen to be 
algebraic and close to 1. The latter is over an energy interval that includes a large number 



of pairs (i, f) and is, in fact, necessary due to large fluctuations in Lif . In other words, 
spectral averaging proves essential both in order to derive a closed form analytical expression 
and to ensure convergence of numerical results. Clearly, an underlying assumption is the 
validity of the "ergodic hypothesis" - that the two averages are equivalent. 



if 



2 

is 



This Section is organized as follows. First, we show that the magnitude of 

determined by the hierarchy of (odd) pairs (k x — n x ) and (k y — n y ), where i = (n x , n y ) and 

^ 2 

/ = (k x ,k y ) are the nearest energy levels. For a given pair, Lij cx H with the largest 
coefficient, by two orders of magnitude, corresponding to k x = n x ±l, k y = n y ^fl. Moreover, 
we show that the probability of the latter oc A/" -1 / 2 , while for most pairs it is oc A/" -1 . 
Consequently, such pairs give an overwhelming contribution to the spectral average, which 
results in eq. (2). This is subsequently verified numerically. 



B. Spectrum and 



in a rectangular box 



For a rectangle, the energy eigenvalues are 

"""»~ 2m Vi| £?/ 
We consider rectangles of the same area A = L x L y in an ensemble with different values of 
the ratio a = L 2 x /L 2 y . Numerically, we use algebraic numbers for a to reduce accidental level 
degeneracies. Expressing the energies in terms of the mean level spacing, 

A = — r 

mA 

we find the dimensionless form of the spectrum 



(4) 



e nxny = \ (a^nl + a^n\) = [n\ + an\) (5) 

For a ~ 1, we have a simplified expression for the spectrum 

e nx n v w \ K + "!) = ^ (6) 

which will be used in the analytical derivation below [9] (in our numerical work, we also use 
algebraic cc's close to 1). For an energy e>l, where the relevant quantities can be described 
semiclassically, the level % nearest to (and below) e will be characterized, in general, by a 
different pair (n x , n y ) for each a. On the average, the level number M of level % = (n x , n y ) is 

(A/") = e > 1 (7) 

In what follows, the variations in Af with a are not important. Consequently, we drop () 
and identify Af with e. 

The matrix element of the orbital angular momentum between the levels i = (n x ,n y ) 
and / = (k x , k y ) (£/>£« for definiteness) is given by 



where (k x — n x ) and (k y — n y ) are odd. Retaining only the terms that contain (k x — n x ) and 
(k y — n y ) , we find a simplified form 

1 





2 16 




7T 4 _ 



-1/2, 



a ' n x a ' n y 



(kx Tlx) \ky T^y) V kx ^ a 

For levels i and / that are nearest in energy 



-, 2 



kl + k 2 y « + nj 



and for small (A^ — ra x ) and (A^ — n y ), we find 

^a; (k x fix) ~ (^j/ ^j/) 



Consequently, for a ~ 1, 



if 



64 



^ (k x -n x ) (ky-riy) 



(9) 



(10) 



(11) 



(12) 



Clearly, the magnitude of 



is determined by the hierarchy of values \k x — n x \, 



"if 



as a functions 



\k y — n y \ such that sgn (k x — n x ) = — sgn (k y — n y ) . Fig. 1 shows 
of M for a single aspect ratio a = 8/31 - 6 « 1.01923. The top straight line corresponds to 
k x — n x = ±1, k y — n y = =)=1 The other two lines correspond, respectively, to k x — n x = ±1, 
k y — n y = =f3 and k x — n x = ±1, k y — n y = =f5 and rr <-> y permutations. The slopes of 
other lines are too small for them to be visible in this plot. Using eqs. (6) and (7), we have 
n x + n y ~ 4A^" Z 71 "- Combining this with (11) and substituting into eq. (12), we find 

for ±l,=Fl 

■4 lor : 1. : 3 and : 1. : 3 (13) 

for ± 1, =f5 and =p 1, ±5 
for the three lines shown in Fig. I. As seen from the figure, these straight lines are in 
excellent agreement with the numerical results. 
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We emphasize that Fig. 1 should be understood as follows: as i moves up, from one 



level to the next, 



Lii 



"jumps" - up or down - between the points on straight lines (whose 
total number is of order A/", with only three shown here), indicating orders of magnitude 
fluctuations as a function of the position in the spectrum. 



C. Spectral average of 
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It follows from (5) that the greater the range of a, the greater is the spectral (x,y) 
asymmetry and the higher in spectrum it is necessary to move to eliminate it 2 . On the 
other hand, keeping a's close to a fixed value introduces a problem of sampling, which, 
again, requires higher energies to be considered numerically. 3 Thus, there exists an inherent 
technical difficulty with ensemble averaging in a box. As a practical matter, we take up to 
400 values of a G [1, 2] for the energy range e ~ 10 4 — 10 7 . This choice of parameters proved 
suitable for ensemble averaging in a numerical part of the study of statistical properties of 
the energy spectrum itself. [8] However, ensemble averaging with these parameters does not 



lead to the numerical convergence of 



L 



if 



A portent of this can be already glimpsed from the orders of magnitude fluctuations in 
discussed in the preceding subsection. Moreover, if M is fixed instead but different 



Li i 



are as large as a function of a as they 



aspect ratios are considered, the fluctuations of 
are as a function of Af. Consequently, together with a-averaging, an additional spectral 
averaging over the energy interval [e — E/2, e + E/2] , E <C e, is performed to achieve 
numerical convergence. The interval width is taken to be > y/e, which, for the above 
parameters, increases the sampling range by several orders of magnitude. Interestingly, it 
also turns out that spectral averaging is amenable to an analytical derivation, which we 
proceed to outline here. 



L if 



is dominated by the contribu- 



The key idea in this derivation is that the average of 
tion from the top line in Fig. 1. Numerically, it is clearly seen from Fig. 2 where the higher 
set of points corresponds to averaging over all pairs (i, f) (all lines in Fig. 1), while the lower 



2 In other words, changes of n x vs. changes of n y lead to disparate changes of energy for rectangles 
whose sides are substantially different. However, this anisotropy is eliminated for an algebraic 
aspect ratio once above a sufficiently high energy in the spectrum. 

3 That is, at low energies the spectra for different o's are very close to each other. 
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set corresponds to averaging only over the top line in Fig. l(k x — n x = ±1, k y — n y = =Fl). 
The fitting lines for the two sets are given, respectively, by 



0.177V A/" 



(14) 



all 



and 



= 0.172V7V 7 



(15) 



top 



which are remarkably close. Note that we used a combined average over 400 values of 
a G [1,2] and over energy window E = 4 x 10 4 . Great sensitivity to the spectral averaging 
is obvious from Fig. 3 where the interval was reduced to E = 4 x 10 2 . 

Physical interpretation of these results is as follows. As has been mentioned in the pre- 
ceding subsection, at a spectral point M there are ~ Af lines whose slopes can be determined 
(for small \k x>y — n Xjy \) via a procedure that resulted in eq. (13). For a given a, as i moves 



if 



"jumps" between these straight 



upward from level to level (that is, as N is increased), 
lines. In principle, a "jump" can be from any one of these lines to any other. However, in 
determining the average, one must remember that not only the slope decreases rapidly from 
the top line down, but the probability of being on a line also decreases rapidly from ~ A/"~ 1//2 
on the top line to ~ A/" -1 over the course of \fN lines. Given this, and in view of the nu- 



merical evidence above, we approximate 
Fig. 1 alone. From (13), we find 

128A/" 



if 



by the contribution from the top line in 



'if 



7T 



128A/" 



P (k x -n x = ±1, k y -n y = =Fl) 
P (k x -n x = ±1) P (k y -n y = ^l\k x -n x = ±1) 



(16) 
(17) 



Here P (k x — n x = ±1, k y — n y = =Fl) and P (k y — n y = =)=1 | k x — n x = ±1) are, respec- 
tively, the probability that the nearest energy level pair satisfies the condition (i, f) = 



(k x — n x = ±1, k y — n y = =)=1) and a conditional probability that k y 



n„ 



=Fl, given that 



k x — n x — ±1. In what follows, these are evaluated in a series of consecutive steps. 
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First, we derive the probability distribution function p (n x ). Consider the energy interval 
defined, in accordance with (6), by the numbers A 2 > Ni >> 1. A narrow interval is defined 
as such that 



AiV AM 

where A is the center of the interval and AN is its width: 

Ao + At 4 4 4 4 

A = 2+1 =-e = -M, AN = N 2 — Ni = —E = -AM 

2 IT IT TT TT 



(18) 



(19) 



All pair points (n x ,n y ) in this interval are located between the two quarter-circles shown in 
Fig. 4, 



Ai <n 2 x + nl < A 2 , n x , y > 



(20) 



We propose a simple ansatz whereby p (n x ) dn x is just the number of states in the interval 
dn x , as illustrated by the shaded areas, which gives 



v M = 



ttAA 



V 



n: 



I 



\ 



0<n x < v^/Vi 
/AT < n x < x/As 



(21) 



Fig. 5 shows this formula vis-a-vis the numerical evaluation for the same a as in Fig. 1 and 
for 0.9 x 10 5 < e < 1.1 x 10 5 . Clearly, the two are in excellent agreement. 

Next, we evaluate the distribution function p (k x — n x ) by convoluting p (n x ) and p (k x ) 
under the assumption of no correlations between n x and k x : 



P (k x -n x )= p (k x -n x + t)p (t) dt 

Jo 



(22) 



The resulting formula is very complicated and in Fig. 6 we will only show its plot vis-a-vis 
numerical evaluation for 0.99 x 10 6 < e < 1.01 x 10 6 for a single a and for an ensemble 
average over 100 algebraic a G [1,1.25]. Clearly the agreement between the two is very 
good. Notice that the dip at zero in the numerical distribution function is due to the fact 
that the probability of k XjV — n Xjy = is suppressed for the nearest levels (see, for instance 
eq. (11) showing that (k x — n x ) and (k y — n y ) should be finite and of opposite signs). Notice 



8 



also that we treat (k x — n x ) as a continuous variable and so a deviation from the numerical 
results should be anticipated for small discrete values of this variable. 

As is clear from (17), we only need the probability P (k x — n x ± 1) for our purposes. 
However, for the reasons just mentioned, we cannot precisely determine it from the distri- 
bution obtained via (22) and shown in Fig. 6 by solid line. Consequently, we use an ansatz 
where this probability is approximated by the area of width ~ 1 near the maximum of the 
distribution, that is 

P(k x -n x = ±l) »p(0) (23) 

where 

32 l + (l-<5) 3/2 + 8K(l-5)-(2-5)E(l-5) 

and K and E are the elliptic functions. Expanding this expression for small 5, while simul- 
taneously replacing N 2 by N, we find 

, x 3 + 81n2-21n5 
p « = 25 



Fig. 7a shows numerical results for P (k x — n x = ±1) vis-a-vis eqs. (24) and (25). Fig. 7b 
also shows P (k x — n x = ±1) \/JJ vis-a-vis p (0) VJ\f. The latter eliminates the main \fM 
dependence and shows that the remainder is a slow-growing function of Af. Unfortunately, 
due to indeterminacy inherently present in our ansatz, it is incapable of exactly describing 
this function. 

To evaluate the conditional probability P (k y — n y = =pl | k x — n x = ±1), we notice that 
by substituting k x — n x = —1 into (10) and using (6) one finds 4 

k y - n y = ^N-(n x -l) 2 - ^fN-nl (26) 

whereof 



4 Evaluating specifically P (k y — n y = 1 | k x — n x = — 1). 
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1 

n x = - 



1 + (ky ~ Uy) 



^N-{k y -n y f 

l + (ky- riyf 



(27) 



in the limit N ^> 1. Consequently, using the distribution function (21), we find the distri- 
bution function for p (k y — n y \ k x — n x = —1) as 

p(k y -n y \k x -n x = -1) =p (n x (k y - n y )) — — - — - (28) 

a [K y — n y ) 

where is n x (k y — n y ) is a function given by eq. (27). Combining now eqs. (26)-(28) with eq. 
((21) we find the analytical form (not shown here due to its complexity) of the conditional 
probability distribution function and plot it in Fig. 8 vis-a-vis the numerical results. 

Before discussing analytical vs. numerical results, we wish to pause briefly on the slight 
difference between p {k y — n y \ k x — n x — — 1) and p (\k y — n y \ \ k x — n x = 1), as seen in Fig. 
8. It follows from a more careful examination of eq. (10). Indeed, it should be written as 

kl + tfnnl + nl + A (29) 

where A > is the separation between the levels. Consequently, (11) should be replaced by 

n x (k x - n x ) pa -n y (k y - n y ) + ^ (30) 

Obviously then, it is possible to have k y — n y = when k x — n x = 1 but not when 
k x — n x = —1. This is reflected in numerical curves in Fig. 8 where the condi- 
tional probability p{k y — n y \ k x — n x = — 1) does not have a value at k y — n y = while 
p (\k y — n y \ \ k x — n x — 1) does. 

We now turn to a closer analysis of Fig. 8. Obviously, the numerical distribution is much 
broader than the analytical distribution. The reason for that is that the former was obtained 
using a large number of a -values whereas the latter was obtained in the approximation where 
a was set to ~ 1 and hence (6). Clearly, even if the a- values are sufficiently close to 1, as in 
our simulations, different ensembles are nonetheless sampled for different a's (for sufficiently 
high energies, as has been previously mentioned) using the exact formula (5). As a result, 
the likelihood of larger \k y — n y \, given that \k x — n x \ — 1, is higher in this case than in our 
analytical derivation based on (6). 

10 



Of note, however, is that for small \k y — n y \ the analytical curve is well approximated by 

v(\k -n\\\k n I - 1) ~ 2 I AN Ifc ^ l= ° 2 I AN (31) 

p(\k y n y \\\k x ^|- 1 )~ 7r(l + fe _ %|2) +2^v ^ n + 2^N (31) 

We assume that, up to a constant, this gives a good approximation to the dependence 
on AN/N = AAf /Af. Consequently, we propose, similarly to p(k x — n x ) before, that 
P {k y — n y = =Fl I k x — n x = ±1) can be approximated by 

P (k y — n y = =Fl I A;^ — ria. = ±1) pa const x p (0 | |fc x — n x | = 1) (32) 

Fig. 9 shows a plot with an empirical const = .225 used both in "exact" p (0 | |fc x — n x | = 1), 
obtained from (28), and approximation (31) vis-a-vis the numerical result. Clearly, it pro- 
vides credence to our ansatz. 



D. Summary 



The key results of this section are summarized in Figs. 1 and 2 and eqs. (14) and (15): 



U 



if 



will fall on one of 



• Given the energy (the level number) ~ Af, the magnitude of 

~ Af lines, such as the three straight line shown in Fig. 1 and given by eqs. (13). 
Which line specifically it will be on depends on a particular aspect ratio of the rectangle 
a. Conversely, for a given a, the line it will be on depends on a specific nearest levels 



pair (i, f) in the vicinity of the energy considered. In other words, 



L 



if 



experiences 



orders of magnitude fluctuation both as a function of a and as a function of energy. 

For numerical convergence it was necessary to perform a combined ensemble aver- 
aging (over a) and energy averaging (over energy interval AAf <^ A/"). The result 
is given by eq. (14) and its magnitude almost entirely derives from the contribu- 



tion of the top straight line in Fig. 1, 



'if 



with (15). The probability to find 

2 

if 



Li 



if 



why 

P(ky 



= 128A/"/7T 5 , as seen by comparison 
on the top line is oc A/" -1 / 2 , which explains 
oc VJ7. It is central to our derivation that the conditional probability 
=Fl I k x - n x = ±1) only weakly depends on A/", ~ 0.14 (1 + AAf /4Af). 
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II. ORBITAL SUSCEPTIBILITY OF A FREE-ELECTRON GAS 



A. Introduction 

It was recently proposed that, at zero temperature the orbital magnetic susceptibility, 
of free electrons in disordered systems can be explained by the two level van Vleck response 
that involves the last occupied (Fermi) level and the first unoccupied level [1], [10]. Whereas 
for an occasional Fermi level in a given realization of disorder, or for an occasional real- 
ization of disorder given a Fermi level close to a fixed energy value, the response can be 
diamagnetic, in the vast majority of cases it is paramagnetic. This prediction was recently 
verified numerically in Ref. [11], which also confirmed that both the mean susceptibility and 
the susceptibility distribution function (mesoscopic fluctuations) can be quite accurately 
described by the two level model [1], [10]. 

Orbital susceptibility of rectangles has been previously studied both analytically and 
numerically [12]- [14]. However, analytically the approach that generally works only at 
sufficiently high temperature was used and numerically only the rational aspect ratios of 
rectangles were examined (e.g. squares). Conversely, here we are concerned with a strictly 
T = response for the irrational (algebraic in this numerical evaluation) aspect ratios. Our 
main conclusion is that, in complete analogy with the disordered systems, both the average 
and the fluctuations are successfully described by the two level van Vleck response that 
involves the Fermi level and the first unoccupied level. Furthermore, the largest contribution 



to the response comes from the values 



Lif 



from the top straight line in Fig. 1, as explained 
in the preceding Section, while the largest fluctuation occur when, for the points on that 
line, the energy difference Sf — Si between the nearest levels is particularly small. 

It will be observed that the orbital susceptibility exhibits a striking absence of self- 
averaging. While in part it is due to the fact that we evaluate the zero-field response at 
zero temperature, the underlying physics underscores, nonetheless, that mesoscopic effects 
are much more pronounced in classically integrable than in classically chaotic systems. 



12 



B. Orbital susceptibility of a rectangular box 



Using the dimensionless notations, where the energy is measured in units of A and the 
susceptibility in units of /^/A, we find the total orbital susceptibility as follows: 



Xtot 



2tt (i\x 2 + y 2 \i) 



* oo 2\L if 

+E E 



i=i 



=1 f=Af+l 



Ef - Si 



(33) 



where £j, % = (n x , n y ), is the unperturbed (zero-field) spectrum (5). While the Landau gauge 
was used in this expression, the final result is gauge-independent (for a discussion, see Ref. 
[11]). The diamagnetic matrix elements are easily calculated and are given by 



(i\x 2 + y 2 \i) 



A 
12 



a 



1/2 1 _ 



nni 



+ a" 1 ' 2 1 - 



Tin: 



(34) 



Fig. 10 shows the result of numerical evaluation of (xtot) plotted as a function of Af, vis- 

2 



a-vis ( 2 



if 



(ef — Ei)' 1 ^, where £j and Ef are now limited to the Fermi level and the first 
unoccupied level, i = Af, f = Af+ 1. This is motivated by the surmise that the contributions 
of the two sums in eq. (33) - diamagnetic and paramagnetic - largely cancel each other over 
the Fermi sea and the total susceptibility, on average, can be explained by a single term 
in the van Vleck sum, namely, the one between the last occupied and the first unoccupied 
levels. (The analogous surmise in disordered systems [1], [10] had been already verified 



1 2 L if (ef-Ei) M , wi 

\ / top 



with the contributions 



numerically [11].) The subset of the latter, ^2 
only from the top line in Fig. 1 is also shown. While the difference in the distribution 
function of (e/ — £*) on the top line relative to the Poissonian should be noted (and will 
be discussed in a separate publication), the dominance of the k x — n x = ±1, k y — n y = =Fl 



contribution to 



Li 



suggests that the top line should dominate the contribution from the 



if 



((£/ — £j) : ), is also shown in Fig. 10, 



% = Af, f = M + 1 term also. Finally, 2 

where / Lif \ is shown in Fig. 2 and approximated by eq. (14) and ((£/ — £j) _1 ) is 
evaluated numerically and found to be (in units of A -1 ) 

1 



£ f 



15.5 



(35) 
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Despite a combined averaging over 400 cc's and energy interval 4000-wide, the absence 
of self averaging is still evident in this plot. (To further emphasize the predominance of 
large fluctuations, in Fig. 11 we show the same (xtot) as in Fig. 10 vis-a-vis (xtot) that 
was obtained for the same a-ensemble but whose energy averaging was performed over 
intervals 10 times narrower.) On the other hand, our surmise that the two level van Vleck 
paramagnetism accurately describes the average orbital response is evident from Fig. 10 as 
the structure of (xtot) vs. Af is well reproduced by the nearest level contributions alone; the 
difference between the two are the contributions from the terms in the double-sum of (33) 
that are due to the levels further below and above the Fermi level. 

The large value of ((£/ — £i) _1 ) in (35) is readily understood from the exponential (Pois- 
sonian) distribution function of the level spacings [5], [6], p(sf — Si) = exp [— (ef — £,)], 
whereof 



where e is a cut-off. Since we evaluate the zero field response at zero temperature, the 
cut-off is the smallest spacing observed for the values of a and energies considered here, 
which happens to be ~ 10~ 7 (the mean level spacing being 1). As already mentioned, 
the distribution function of the level spacings on the top line of Fig. 1 will be discussed 



out however that in reality, even at zero temperature, the magnetic field itself introduces a 
natural cut-off; for disordered systems this was discussed in Ref. [10]. 

We now turn to the study of fluctuations. To understand the nature of mesoscopic 
fluctuations, we investigate the energy differences between the nearest levels. We limit 
our consideration to the top line in Fig. 1 as it dominates the response, both in terms 
of the mean and the fluctuations. Figs. 12-14 show plots of (ef — and 



function of M for three particular values of a, respectively a,\ = 457 a855 /185 ~ 1.01638, 
a 2 = 911 a755 /150 « 1.14379, and a 3 = 643 a655 /63 « 1.09657. First, we notice the existence 
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(36) 



elsewhere, however, it turns out that ((£/ — £j) 1 ) is also formally divergent and, for the 
parameters used here, numerical evaluation gives ((£/ — £i) _1 ) ~ 12.5. It should be pointed 



of very large values of (e/ — e^ 1 , which are the source of very large mesoscopic fluctuations. 
Second, we notice patterns in the structure of (sf — Ej) 1 and Sf — E{ as a function of Af. 
Lastly, we notice that the size of the peaks and the pattern structure is very strongly in- 
dependent. 

It turns out that the patterns in Figs. 12-14, as well as the size and the location of the 

peaks, can be determined analytically. Considering for simplicity the case of a — 1 + f3, 

< < 1/2, we find 

4« 1/2 _ f2(n x - n y + 1) - (3(2n y - 1)\ (k x - n x = 1, k y - n y = -1\ 

vr [£f £i) ~ \2(n y -n x + l)+(3(2n y + l))\k x -n x = -l,k y -n y = lJ 

The algebraic (5 can be represented as a series 

111 

13 = - + - + - + ... (38) 
p q r 

where p > 2, \q\ > p 2 , \r\ > q 2 ,. . . are integers. The series can be truncated 
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1 








+ v 
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P = 


1 
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1 




+ - 


+ - + V 


P 


Q 


r 



(39) 



where rj is the residual algebraic number. For /3's corresponding to Figs. 12-14, the truncated 
series are as follows: 

A = 5.1 x 10" 11 

1 61 81672 

^-^l4 + 2^460 + 5 - 6 x 10 " 13 ^ 

/3 3 = — - l — x 10- 13 

1 10 292 472427 

Specifics of truncation depend on the position in the energy spectrum, to which numerical 
analysis is extended, and on the values of integers p, q, r, . . ., as explained below. 

The key to understanding the spacings structure in Figs. 12-14 is that it can be com- 
pletely described by as few as the first two rational numbers in the approximation of j3 (39). 
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(Only for the very rare occurrences of degeneracy, Sf — Si = 0, does the residual term r\ needs 
to be considered.) Therefore, we turn to the analysis of the interplay between the rational 
numbers. We begin with the spacings structure for (5 = p -1 , shown in Figs. 15 for p = 100. 
Analytically, (37) yields the following series of straight lines, as a function of n y : 

(ef-ei) =2C + 2+ ny, r " o 2 ),k x -n x = l,ky-n y 



7T v 1 l/ p p y ' V C = 0,1,2,... 

(ef-Ei) =-2C + 2 + - + -n y l yK ' 2 ~ 2 ),k x -n x = -l,k y 



P P V C = 1,2,3, 

Here C = n x — n y is a non-negative integer, due to the choice of a > 1, and the limits on n y 
are determined from the condition 

4a 1 / 2 

0<^— (e/-£i)< 2 (42) 

7T 

where the second inequality can be understood from the fact that changing C by 1 changes 
4a 1 / 2 (ef — Si) /it by 2, making it extremely unlikely that Ef is the nearest level to £j if the 
second constraint in (42) is not satisfied; this is also confirmed numerically. Furthermore, 
the maxima in Fig. 15a (minima in Fig. 15b) are found as follows: 

( £ f ~ £ mi „ = when ( n _ n 1 o ),k x -n x = l,k y -n y = -l 

(43) 



4 « V % ^ 1 i, fn y =p(C-l)\ . 1 . 

(£/ — £j) = -, when " , K;,. — n x — — 1, k v — n v — 1 

n K 1 lJ p \C = 2,3,4,.../' '2/2/ 

Combining eqs. (41) and (43) with 

A^^Sl ,44, 

completely describes the spacings structure as a function of A/", for instance, as shown in 
Figs. 15. 

We turn to the next approximation in (39), f3 = p~ 1 + q~ 1 . For simplicity, we will consider 
only k x — n x = 1, k y — n y — —1 and q < 0. In this case, the spacings structure is described 
by the following equations: 
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4a 1 / 2 2t + 2p+l l-2n y ( 1 < % < - p < t < \ 

vr ^ £iJ p q ' < % < (g^te t < -p - J 1 J 

Here t is an integer and possible n y 's are subject to the constraint 

t = pC-n y (46) 
where, as before, C = n x — n y > is an integer. The minima of the structure are found at 

-t,t<-p-l (47) 



n y =p 



t | (2t + 2p+l)g h J_ + 1 
p 2p 2 2p 



where the [J brackets denote the floor (integer value) function. Together, eqs. (44)- (47) 
completely describe the spacings structure in Figs. 12. Furthermore, generalization to the 
next iteration, (3 = p L + q^ 1 + r _1 , is straightforward and accounts for Figs. 13 and 14. 

Semiquantitatively, the peak structure of (sf — e^)" 1 along Af axis is summarized below. 
In what follows, we will discuss the spacings structure as a function of 2n y , which is easily 
converted to that of A/", as explained above. The upper limit of Af, vs. the values of p, q, 
r, ... in approximation of a particular a ((3), determines how many terms is necessary to 
keep in (39). Namely, q~ x enters in (39) when 2n y > q/p and r -1 enters when 2n y > r/q; 
this is illustrated by (40), as applied to Figs. 12-14. Noting that peaks correspond to the 
end points (negative slope) and first points (positive slope) of the straight lines given by 
eqs. (41), (45), etc., and calling the distance between the tallest peaks "period," we find the 
following 5 : 

• Peaks of p _1 : "period" ~ p; height of peaks is p. 

• Peaks of p -1 + q~ l : (assuming that pq and p + q have no common factor) period ~ 2pg, 
with ~ p 2 peaks per period separated by distance ~ 2g/p; upper limit of peaks is pq 



5 It should be emphasized that not all peaks are observed for the Af considered; in fact it is possible 
to not observe any peaks at all, in which case the maxima of the structure are determined by the 
end points of the interval on the Af axis. 
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(heights of peaks are pq, pq/3, . . ., ~ q/p, or in reverse order), for p + q odd, and oo 
(oo, pq/2, pq/4:, . . . , ~ q/p) for p + q even (oo indicates level degeneracy, s$ — £j = 0, 
so that the residual term r\ must be considered in (39)). 

• Peaks of p~ l + g -1 + r _1 : (assuming that pgr and pg + gr + pr have no common factor) 
period ~ 2pqr, with ~ p 2 q 2 peaks per period separated by distance ~ 2r/pq; upper 
limit of peaks is pqr (heights of peaks pqr, pgr/3, . . ., ~ r/qp, or in reverse order), for 
pq + qr + pr odd, and oo (oo, pqr/2, pqr/4, . . . , ~ r/qp) for pq + qr + pr even. 



C. Summary 

The central result of this section is shown in Fig. 10. It demonstrates that the mean 
zero-temperature, zero-field orbital magnetic susceptibility of a free electron gas in a rect- 
angular box can be explained in terms of a two-level van Vleck response - that of the last 
occupied (Fermi) and the first unoccupied levels. Furthermore, it is dominated by the con- 
tributions from the top line in Fig. 1, namely k x — n x = ±1, k y — n y = =1=1, which is also the 



largest contributor to 



L 



if 



y '"y 

2\ 



discussed in the preceding section. In fact, the mean value 



,2 

L if 



(( £ f~ £ i) X >, where ((£/-£*) 



of susceptibility is reasonably well described by 
is large due to the absence of correlations for small level separations. 

It is also evident that the orbital susceptibility is largely a non self-averaging quantity, 
as seen from Figs. 10 and 11. This is due to the existence of huge variations in inverse level 
spacings, which, in turn, allow for such large contributions that may singularly outweigh 
the totality of more typical contributions in the average response. Such variations were 
explained in terms of a decomposition of algebraic aspect ratios into rational numbers, whose 
interplay in (39) is crucial for understanding the peaks of (Ef — £j) -1 . It must be borne in 
mind, however, that this feature of the orbital susceptibility is very fragile with respect to 
perturbations and that mesoscopic fluctuations will be suppressed at finite temperatures (or 
even by finite values of the magnetic field); we intend to address this problem elsewhere. 



18 



III. CONCLUSIONS 



The most striking feature revealed here is the non-self-averaging property of physical 
quantities in a rectangular box, which represents a class of integrable billiard problems. We 
had previously discussed [8] the persisting long-range correlations in the semiclassical energy 
spectrum of this system. These correlations are more complex than those in classically 
chaotic (disordered) systems [5]- [7]. In particular, we discussed the large, non-decaying 
oscillations of the level number variance on an energy interval as a function of the interval 
width. Similarly, we find that mesoscopic fluctuations here are much more pronounced than 
in classically chaotic systems. For instance, while eqs. (1) and (2) point to the same order 
of magnitude in a ballistic disordered system, t ~ L, and a rectangular billiard, the latter 
will have much larger fluctuations (we have discussed the difficulties with averaging in text). 

The one similarity that stands out for both integrable and chaotic systems is that both 
the average orbital susceptibility of the free electron gas and its fluctuations can be well 
described by a two level van Vleck response that couples the last occupied (Fermi) and the 
first unoccupied levels. For disordered systems, this has been demonstrated previously in 
Refs. [1], [10], and [11] and for an integrable case in this work. The difference, however, is 
that in disordered systems the non-self-averaging effects are less pronounced: in the absence 
of cut-offs (temperature, finite magnetic field, etc.), the average is well defined and only the 
higher cumulants are divergent. In a rectangular box, even the average is already ill-defined, 
as pointed out in discussion of eq. (36). 

Our next step will be to investigate the effect of temperature on the orbital magnetism 
of integrable systems. Towards this end, we will apply Imry's formalism, which allows to 
express the average response in terms of the level correlation function. [15] For rectangles, 
the latter is now well understood, including as function of magnetic field. [8] This formalism 
works well at sufficiently high temperature and should provide an insight into the scales at 
which the transition to the zero-temperature limit, discussed here, occurs (as was done for 
disordered systems [1]). Our results will then be discussed vis-a-vis previous works [13], [14]. 
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IV. FIGURE CAPTIONS 



1. Figure 1 



of 

0, 



Lif 
Lif 



L; 



if 



vs. Af (e — Af in our approximation). The number of lines, where the value 
might fall, increases with Af '. For a given Af (more precisely, for a given level 

2 

is on one of the lines, depending on a. As, for any given a, Af is increased 



(i — > i + 1 = /, / — > / + 1), Lj/ "jumps" to another line. The lines are the same regardless 
of a. For illustration, we present numerical results for a = 8/31 ' 6 ~ 1.01923. Equations for 
the three straight lines shown here are given by (13). 



2. Figure 2 



Numerical average of 



Lif \ vs. Af vis-a-vis the fits (14) and (15). A combined average 
over 400 algebraic values of a G [1,2] and over the energy interval [Af — 2 x 10 4 , Af + 2 x 10 4 ] 
was used. The bottom line and fit correspond to the contribution from the top line only in 
Fig. 1 (k x -n x = ±1, k y -riy = =Fl). 



3. Figure 3 

Same as the top set of dots in Fig 2 vs. the result of averaging over a narrower interval, 
[Af -2 x 10 2 ,AT + 2 x 10 2 ]. 

4- Figure 4 

Shaded areas represent the probabilities p (n x ) dn x . 

5. Figure 5 

Distribution function p (n x ) vs. n x : analytical result with p (n x ) given by eq. (21) 
vis-a-vis numerical data for 0.99 x 10 5 < e < 1.01 x 10 5 and the same a as in Fig. I. 
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6. Figure 6 



p(k x — n x ) vs. (k x — n x ): analytical result obtained via (22) vis-a-vis numerical data for 
0.99 x 10 6 < e < 1.01 x 10 6 for the same a as in Fig. 1 and for an ensemble average over 
100 algebraic a G [1,1.25]. 

7. Figure 7a 

Probability P (k x — n x = ±1) vs. M: analytical ansatz based on (23) and (24) and on 
its approximation (25) (solid lines) vis-a-vis numerical data. 

8. Figure 7b 

P (k x — n x = ±1) VJ7 vs. Af: same as Fig. 7a, but with the main dependence, oc A/" -1 / 2 , 
eliminated. 

9. Figure 8 

Distribution functions of conditional probability p (k y — n y \ k x — n x — — 1) and 
p(\k y — n y \ | k x — n x = 1): analytical result obtained via (28) vs. numerical data. Notice 
that for the latter p (0 | k x — n x = —1) = while p (0 | k x — n x — 1) ^ 0, as explained in 
text. 

10. Figure 9 

Conditional probability P (k y — n y = =Fl | k x — n x — ±1) vs. M: an ansatz based on 
0.225 x p (0 | \k x — n x \ = 1) and its approximation 0.14 x (1 + AN/4N) (solid lines) vis-a- 
vis numerical data. 
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11. Figure 10 



Total magnetic susceptibility, evaluated via (33) and (34) - top set of dots. 



>-i 



where £j and £/ are the Fermi level and the first unoccupied level 



respectively - middle set of dots. ( 2 



L 



if 



(sf — Si) 1 ) , where £j and £/ are the Fermi 

/ top 



level and the first unoccupied level such that k x — n x = ±1, k y — n y = =Fl (top line in Fig. 1) 



bottom set of dots. 2 



given by the top set of dots in 



((£/ - Si)' 1 ) with ^ L, 
Fig. 2 and ((£/ — £i) _1 ) ~ 15.5. A combined average over 400 algebraic values of a e [1,2] 
and over the energy interval [A/" - 2 x 10 3 , A/" + 2 x 10 3 ] was used. 



12. Figure 11 

Large dots same as top dots in Fig. 10. Small dots - a combined average over same 400 al- 
gebraic values of a G [1, 2] but over the narrower energy interval [Af — 2 x 10 2 , Af + 2 x 10 2 ]. 



13. Figure 12a 

(s f - Si)' 1 vs. Af for a x = 457 a855 /185 w 1.01638, where e { and £ / are the nearest levels 
such that k x — n x — ±1, k y — n y — =Fl. /5 = q; — 1 is decomposed according to the first 
eq. (40). Coordinates of point A are (718964, 25443.6) and can be described analytically by 
(45) with t = -62, C = n x - n y = 10 and n y = 672. 



1^. Figure 12b 

Same as Figure 12a, shows Ef — £j vs. A/" and point A. 



15. Figure 13a 

(sf-SiY 1 vs. A" for a 2 = 911°- 755 /150 w 1.14379 and /3 given by second eq. (45) 
Points A-D have been examined; they can be described analytically by equations that are 
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an extension of (45) to include the next iteration in (39). 



16. Figure 13b 

Same as Figure 13a, shows £/ - £j vs. Af and points A-D. 

17. Figure l^a 

(e f - Si)' 1 vs. Af for a 3 = 643 a655 /63 « 1.09657 and (3 given by third eq. (45). Points 
A and B have been examined; they can be described analytically by equations that are an 
extension of (45) to include the next iteration in (39). 

18. Figure 14b 

Same as Figure 14a, shows £/ - e ; vs. Af and points A and B. 

19. Figure 15a 

(ef — Si) 1 vs. Af for the rational aspect ratio a = 1.01. All the points are described 
analytically by eqs. (41) with p = 100. The maxima, denoted by A and B (two points each, 
corresponding to increasing and decreasing functions of Af in (41)) are given by eqs. (43). 

20. Figure 15b 

Same as Figure 15a, shows ej-£j vs. Af and points A and B. 
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